1. Field of Invention
The present invention relates to computerized simulation of hydrocarbon reservoirs in the earth, and in particular to simulation of fluid movement and other state changes in hydrocarbon reservoirs formed of layered formations.
2. Description of the Prior Art
In the oil and gas industries, the development of underground hydrocarbon reservoirs often requires the building of computer simulation models. These underground hydrocarbon reservoirs are often complex rock formations which contain both a petroleum fluid mixture and water, with the reservoir fluid content, often existing in two or more fluid phases. The petroleum mixture is produced by wells drilled into and completed in these rock formations.
Reservoir simulation belongs to the general domain of simulation of flow in porous media. However, reservoir simulation normally involves multiple hydrocarbon components and multiple fluid phases in an underground geological formation which is under high pressure and temperature. The chemical phase behavior of these hydrocarbon fluids and the included groundwater has to be taken into account in these simulators.
Sometimes, fluids such as water and/or gases are also injected into these rock formations to improve the recovery of the petroleum fluids. Simulation models therefore contain data which describe several types of information: the specific geometry of the rock formations and the wells, the fluid and rock property data, as well as production and injection history pertaining to the specific reservoirs of the oil or gas field in question.
Oil and gas companies have come to depend on reservoir simulation as an important tool to enhance the ability to exploit a petroleum reserve. Simulation models of reservoirs and oil/gas fields have become increasingly large and complex. The simulator (known in the petroleum industry as a reservoir simulator) which in the past has run these models was a computer operating under the control of a set of computer instructions or software. The software was in coded form, including some specific numerical algorithms and data constructs of an underlying mathematical model. The mathematical model which represented the physics of fluid movements in these hydrocarbon reservoirs was a system of nonlinear partial differential equations which described the transient multiple-phase, multiple-component fluid flow and material balance behaviors in these reservoirs. The fluid flow and material balance changes were induced by the production and/or injection of fluids as well as the pressure-volume-temperature (PVT) relationships of the reservoir fluids.
The reservoir simulator simulated the material balance and fluid flow in underground reservoirs and the included surrounding porous rock formations by subdividing the volume into contiguous cells known as grid blocks. A grid block was the basic finite volume where the mathematical model was applied. The number of grid blocks needed varied depending on the resolution needed for the simulation and the size of the reservoirs and oil and gas fields in question.
In the structured grid case, the finite difference stencil was typically 5-point or 9-point, in two dimensions, and 7-point or 11-point in three dimensions. The most common case for a structured grid in a three dimensional field-scale simulation was the 7-point stencil. In the semi-unstructured case, the in-layer part of the stencil could have an arbitrary number of connected points, whereas the vertical number of connected points was three. For each time step, the multiple-phase, multiple-component system of the nonlinear discretized material balance equations was typically iteratively solved using what was known as the generalized Newton's method. In the industry, this method has been usually referred to as the Newtonian iterations. At each Newtonian iteration, a linear system of equations was constructed where the matrix, known as the Jacobian matrix, and the right-hand-side vector, known as the residuals, were used to solve for the change in the primary variables of the system.
In the time discretization, when all the primary variables were taken at the new time level of the time step, the method was considered fully implicit (FI). When only the pressure variable was taken at the new time level while all other variables, such as concentration or saturation, were taken at old time level of the time step, it was known as implicit pressure explicit saturation (IMPES). IMPES solved one equation (the pressure equation) per grid block per time step implicitly and was far cheaper computationally than FI, but had stability limitations which constrained the time step size. There were other schemes which adjusted the implicitness dynamically on a cell-by-cell basis, which were known generally as an adaptive implicit method (AIM).
An industry practice for solving this linear system of equations was via a preconditioned iterative method, as the system was normally too big to be solved by a direct method such as Gaussian elimination. The typical state-of-the-art preconditioner used in the industry was nested factorization (NF) or incomplete LU factorization (ILUF) which generated approximate upper (U) and lower (L) triangular matrices. The procedure was highly recursive in nature and was not conducive to parallel implementation.
At the same time, computer hardware has evolved rapidly to become inexpensive and fast. With the availability of fast and inexpensive, commodity-based multi-processor machines such as PC clusters, there was an increasing demand for effective use of these inexpensive systems.
The measure of efficiency of a computational algorithm in parallel computing is its scalability. A method is considered to be perfectly scalable or to have a one hundred percent parallel efficiency if it takes one hour to solve the problem on the computer with a single CPU; 0.5 hours if the work is exactly divided out and given to two CPU'S; the time to solve the same problem using four CPU's is 0.25 hours; and so on. That is, there is no parallelization overhead when a perfectly scalable method is used. This is an ideal situation. In reality, many reasons can cause the solution time on a real system to be far from this ideal.
When the problem was too large to be solved using a single CPU, the grid was decomposed into smaller blocks or chunks, called sub-domains, each to be worked on by a separate CPU. However, this partitioning was artificial and the sub-domains of the grid were in fact connected and communicated with each other. In a reservoir simulation, there must be continuity of flow, pressure, and temperature field across the sub-domain boundaries. Since the native scalar method was recursive and was not easily parallelizable, a multi-level method, for example the additive Schwarz method, or others, was typically used.
The recursive algorithm, like the ILU factorization method, was applied in the local sub-domain and a global iteration step was applied to resolve the discrepancy along the boundaries of the sub-domains. This work to resolve the domain boundaries typically increased as the number of processors (and thus the number of sub-domains) used to solve the system increased. This posed limits on the scalability inherent to this type of extension. Additionally, the convergence of the solution also depended on the additional level of iterations and the number of sub-domains or CPU's used to solve the problem. Hence, the results usually were somewhat different when a different number of CPU's was used to solve the same problem while everything else remained equal. The simulation results for the same model but using a different number of processors were not the same. This could lead to interpretation problems, verification problems, and added uncertainty in the results.
The primary methodology for solving the linear system of the class of simulation problem in question involved nested factorization or incomplete LU factorization, which were recursive algorithms. As discussed earlier, the predominant method for parallel computation was in the form of domain decomposition where a global iterative step was constructed to solve the partitioned system which by the physical nature of the reservoir was tightly coupled. This global iteration was therefore, so far as is known, the standard technique added in the parallel computing environment involving the concurrent use of multiple processors. It was not needed in the serial (single CPU) application. In general, the existing methods had significant dependency on the number of processors used in parallel to solve a problem.
The prior methods used different iterative sequences in parallel simulation as compared to serial simulation. For the same problem, a parallel simulation had an added level of global iterations in parallel processing which increased the processing complexity of the work. This lead to a loss of efficiency as processor count increased.
The simulation results for the same model but using a different number of processors were not the same. This could lead to interpretation problems, verification problems, and added uncertainty in the results.